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Abstract — Attitude control systems naturally evolve on non- 
linear configuration spaces, such as and SO (3). The nontriv- 
ial topological properties of these configuration spaces result in 
interesting and complicated nonlinear dynamics when studying 
the corresponding closed loop attitude control systems. In this 
paper, we review some global analysis and simulation techniques 
that allow us to describe the global nonlinear stable manifolds 
of the hyperbolic equilibria of these closed loop systems. A 
deeper understanding of these invariant manifold structures 
are critical to understanding the global stabilization properties 
of closed loop control systems on nonlinear spaces, and these 
global analysis techniques are applicable to a broad range of 
problems on nonlinear configuration manifolds. 

I. Introduction 

Global nonlinear dynamics of various classes of closed 
loop attitude control systems have been studied in recent 
years. An overview of results on attitude control of a rotating 
rigid body is given in [1]. Closely related results on attitude 
control of a spherical pendulum (with attitude an element of 
the two-sphere S^) and of a 3D pendulum (with attitude an 
element of the special orthogonal group SO (3)) are given in 
[2], [3], [4]. These and other similar publications address the 
global closed dynamics of smooth vector fields. Assuming 
that the closed loop vector field has an asymptotically stable 
equilibrium, as desired in attitude stabilization problems, 
additional hyperbolic equilibria necessarily exist. The do- 
main of attraction of the asymptotically stable equilibrium 
is contained in the complement of the union of the stable 
manifolds of the hyperbolic equilibria. These geometric 
factors motivate the current paper, in which new analytical 
and computational results on the stable manifolds of the 
hyperbolic equilibria are obtained. 

To make the development concrete, the presentation is 
built around two specific closed loop vector fields: one for 
the attitude dynamics of a spherical pendulum and one for 
the attitude dynamics of a 3D pendulum. In analyzing these 
two cases, we introduce new analytical and computational 
tools that are broadly applicable to studying the geometry of 
more general attitude control systems. 

Taeyoung Lee, Mechanical and Aerospace Engineering, Florida Institute 
of Technology, Melbourne, FL 39201 taeyoung@fit.edu 

Melvin Leok, Mathematics, University of California at San Diego, La 
JoUa, CA 92093 mleok@math.ucsd.edu 

N. Harris McClamroch, Aerospace Engineering, University of Michigan, 
Ann Arbor, MI 48109 nhm@umich.edu 

*This research has been supported in part by NSF under grants CMMI- 
1029551. 

tThis research has been supported in part by NSF under grants DMS- 
0726263, DMS-1001521, DMS-1010687, and CMMI-1029445. 



II. Spherical Pendulum 

A spherical pendulum is composed of a mass m connected 
to a frictionless pivot by a massless link of length /. It is acts 
under uniform gravity, and it is subject to a control moment 
u. The configuration of a spherical pendulum is described by 
a unit- vector g G M^, representing the direction of the link 
with respect to a reference frame. 

Therefore, the configuration space is the two-sphere = 
{q G M^|g-g = 1}. The tangent space of the two- sphere 
at g, namely T^S^, is the two-dimensional plane tangent to 
the unit sphere at g, and it is identified with TgS^ :^ {cj G 
I g • = 0}, using the following kinematics equation: 

g = X g, 

where the vector cj G represents the angular velocity of 
the link. The equation of motion is given by 

where the constant g is the gravitational acceleration, and 
the vector 63 = [0, 0, 1] G denotes the unit vector along 
the direction of gravity. The control moment at the pivot is 
denoted by G M^. 

A. Control System 

Several proportional-derivative (PD) type control systems 
have been developed on in a coordinate-free fashion [5], 
[6]. Here, we summarize a control system that stabilizes a 
spherical pendulum to a fixed desired direction g^^ G S^. 

Consider an error function on S^, representing the pro- 
jected distance from the direction q to the desired direction 
qa, given by 

"^{Q^Qd) = l-q-qd' 

The derivative of ^ with respect to q along the direction 
= C ^ where ^ G and ^ • g = 0, is given by 

Dg^(g, gd) • Sq = -(^ xq)'qd = {qd x g) • 

For positive constants kq^k^j, the control input is chosen as: 

u = mf{-k^uj - kqqd x g - yg x 63). 

The corresponding closed loop dynamics are given by 

ib = -kujUJ - kqqd X g, (1) 
q = u) X q. (2) 

This yields two equilibrium solutions: (i) the desired 
equilibrium {q,uj) = (gd,0); (ii) additionally, there exists 



another equilibrium {—Qd^O) at the antipodal point on the 
two- sphere. 

It can be shown that the desired equilibrium is asymptot- 
ically stable by using the following Lyapunov function: 

In this paper, we analyze the local stability of each equi- 
librium by linearizing the closed loop dynamics to study 
the equilibrium structures more explicitly. In particular, we 
develop a coordinate-free form of the linearized dynamics of 
([T]), ([2]), in the following section. 

B. Linearization 

A variation of a curve q{t) on is a family of curves q^{t) 
parameterized by e G M, satisfying several properties [6]. It 
cannot be simply written as q^{t) = q{t) -\- eSq{t) for Sq{t) 
in R^, since in general, this does not guarantee that q^{t) lies 
in S^. In [7], an expression for a variation on is given in 
terms of the exponential map as follows: 



q'{t)=eM4{t))q{t), 



(3) 



for a curve ^{t) in satisfying ^{t) ■ q{t) — for all t. 
The hat map • : R^ ^ 50 (3) is defined by the condition that 
xy = X X y for any x^y G R^. The resulting infinitesimal 
variation is given by 

d 



Sq{t) 



de 



q'{t)=mxQ{t)' 



(4) 



The variation of the angular velocity can be written as 

uj^{t) =uj{t)^eSuj{t), (5) 

for a curve Sw{t) in R^ satisfying q{t) ■ w{t) = for all 
t. Hereafter, we do not write the dependency on time t 
explicitly. 

The time-derivative of dq can be obtained either from ^ 
or by substituting ([3]), ^ into ([2]), and considering the first 
order terms of e. In either case, we have 

Sq = ixq^(^x{ujxq)=Sujxq^ujx{(^xq). 

Using the vector cross product identity a x {b x c) = {a ■ 
c)b — {a ■ b)c for any a^b^c e R^, this can be written as 

ixq^{^-q)w-{^-uj)q = Sujxq^{uj-q)^-{uj- ^)q. 

Since - q = 0, u - q = 0, this reduces to 

^ X q = Suj X q. 

Since both sides of the above equation are perpendicular to 
q, this is equivalent to q x x q) = q x {du x q), which 
yields 

i - {q ■ i)q = q X {^^ X q)- 

Since ^ • = 0, we have ^ • ^ + ^ • ^ = 0. Using this, the 
above equation can be rewritten as 



i = -{C'{uJX q))q^qx {5uj x q) 
= {qq^u)^^{I-qq^)Suj. 



This corresponds to the linearized equation of motion for (|2]). 
Similarly, by substituting (|4]), ^ into ([T]), we obtain 



Suj ■■ 



- kqqd x{^xq) 
kqqdqS,, 

which is the linearized equation for ([T]). 

Equations ([6]), ^ can be written in a matrix form as 

kqqdq 



'i' 




Suj 





-qq 









5uj 



Ax, 



(7) 



(8) 



where the state vector of the linearized controlled system is 
X = 5uj\ G R^. A spherical pendulum has two degrees of 
freedom, but this linearized equation of motion evolves in 
R^ instead of R^. Since q-uo = {) and g • ^ = 0, we have the 
following two additional constraints on ^, 5a;: 



Cx 











"O" 


—uo^ q q^ 




duo 








(9) 



(6) 



Therefore, the state vector x should lie in the null space of the 
matrix C G R^^^. However, this is not an extra constraint 
that should be imposed when solving ([8]). As long as the 
initial condition x(0) satisfies ([9]), the structure of ( p]), ([2]), 
and ([s]), guarantees that the state vector x{t) satisfies^9| for 
all U i.e. j^C(t)x(t) = for alU > when C{0)x{0) = 
0. This means that the null space of C is a flow-invariant 
subspace. 

C. Equilibrium Solutions 

We choose the desired direction as qd = 63. The equi- 
librium solution (qdjO) = (63,0) is referred to as the 
hanging equilibrium, and the additional equilibrium solution 
{—qd^^) = (—63,0) is referred to as the inverted equi- 
librium. We study the eigen- structure of each equilibrium 
using the linearized equation ([8]). To illustrate the ideas, the 
controller gains are selected sls kg = k^j = I. 

1) Hanging Equilibrium: The eigenvalues A^, and the 
eigenvectors Vi of the matrix A at the hanging equilibrium 
(63, 0) are given by 

Ai,2 = (-1 ± V3i)/2, A3,4 = Ai,2, A5 = 0, As = -1, 
vi,2 = ei + (-1 ± ^3^)64/2, ^3,4 = 62 + (-1 ± V^i)e5/2, 

V5 63, ^6 66, 

where G R^ denotes the unit-vector whose i-th element 
is one, and other elements are zeros. Note that there are 
repeated eigenvalues, but we obtain six linearly independent 
eigenvectors, i.e., the geometric multiplicities are equal to 
the algebraic multiplicities. 

The basis of the null space of the matrix C, namely M{C) 
is {ei , 62, 64, 65}. The solution of the linearized equation can 
be written as x(t) = X]i=i ex.p{Xit)vi for constants q that 
are determined by the initial condition: x(0) = Yl^=i^i^i' 
But, the eigenvectors ^5,^6 do not satisfy the constraint 
given by ([9]), since they do not lie in J\f{C). Therefore, 
the constants c^^cq are zero for initial conditions that are 
compatible with We have Re[Ai] < for 1 < z < 4. 
Therefore, the equilibrium {q^uj) = (63, 0) is asymptotically 
stable. 



2) Inverted Equilibrium: The eigenvalues A^, and the 
eigenvectors Vi of the matrix A at the inverted equiHbrium 
(—63,0) are given by 

Ai,2 = -(V^ + l)/2, A3,4 = (V^ - l)/2, A5 = 0, As = -1, 
= ei - (V^ + 1)64/2, ^2 = 62 - (V^ + 1)65/2, (10) 
^3 = (V5 + l)ei/2 + 64, ^4 = (V^ + 1)62/2 + 65, 

^5 = 63, ^6 = 66- 

The basis of M{C) is {61,62,64,65}. Hence, the eigenvec- 
tors ^5,^6 do not lie in J\f{C). Therefore, the solution can 
be written as x{t) = Xlt=i ^« exp(Ait)vi for constants ci that 
are determined by the initial condition. 

We have Re[Ai,2] < 0, and Re[A3,4] > 0. Therefore, 
the inverted equilibrium {q^uj) = (—63,0) is a hyperbolic 
equilibrium, and in particular, a saddle point. 

D. Stable Manifold for the Inverted Equilibrium 

I) Stable Manifold: The saddle point (—63, 0) has a stable 
manifold W^, which is defined to be 



W'{-es,0) = {{q,uj) G TS^ | lim T\q,uj) 

t^oo 



(-63,0)}, 



where : {q{0)^uj{0)) {q{t)^uj{t)) denotes the flow map 
along the solution of ([T]), The existence of W^{—es^O) 
has nontrivial effects on the overall dynamics of the con- 
trolled system. Trajectories in 63,0) converge to the 
antipodal point of the desired equilibrium (63, 0), and it takes 
a long time period for any trajectory near 63,0) to 
asymptotically converge to the desired equilibrium (63,0). 

According to the stable and unstable manifold theorem [8], 
a local stable manifold 63,0) exists in the neighbor- 

hood of (—63,0), and it is tangent to the stable eigenspace 
£^^(—63,0) spanned by the eigenvectors vi and V2 of the 
stable eigenvalues A 1,2. The (global) stable manifold can be 
written as 

W%-es,0) = U-^"*(^ioc(-e3,0)), (11) 

which states that the stable manifold can be obtained by 
globalizing the local stable manifold Wf^^ by the backward 
flow map. 

This yields a method to compute iy^(— 63,0) [9]. We 
choose a small ball Bs c Wi^^^{—es^O) with a radius S 
around (—63,0), and we grow the manifold VK^(— 63,0) by 
evolving Bs under the flow J^~^. More explicitly, the stable 
manifold can be parameterized by t as follows: 



W%-es^O) = {T-\Bs)} 



t>0' 



(12) 



We construct a ball in the stable eigenspace of (—63,0) 
with sufficiently small radius 6, i.e. Bs c £^f^^(— 63, 0). 



From the stable eigenvectors ^1,^2 at ( 10), £^f^^(— 63, 0) can 
be written as 

^foc(-e3,0) = {{q,oj) eJS^\q = exp(aiei +a2e2)(-63), 

oj = -f{-{V5 + l)/2)(ai6i + 0^262) for ai, 0^2 G M}, 

(13) 



where —q'^ in the expression for u corresponds to the 
orthogonal projection onto the plane normal to as required 
due to the constraint q - u = 0. 

We define a distance on TS^ as follows: 

C^TS2((gi,CJi), ((72,^^2)) = \/^((?l,^2) + W^l -^2||. (14) 

For J > 0, the subset Bs of £^f^^(— 63,0) is parameterized 
by eS^ as 

Bs = {{q,uj) e TS^ \q = exp(aiei +a2e2)(-63), 
oj = -g^(-(V5 + l)/2)(ai6i + 0^262), where 
S 



cos 0, 



l/V^+(V^+l)/2 

a2 = ^ ^ sin i9, for e SH. (15) 

l/V^+(V^+l)/2 ^ 

The given choice of the constants ai, 0^2 guarantees that any 
point in Bs has a distance S to (—63,0) according to the 
distance metric ([14]). 

2) Variational Integrators: The parameterization of the 



stable manifold Ws in (12) requires the computation of 
the backward flow map However, general purpose 

numerical integrators may not preserve the structure of 
the two-sphere or the underlying dynamic characteristics, 
such as energy dissipation rate, accurately, and they may 
yield qualitatively incorrect numerical results in simulating 
a complex trajectory over a long-time period [10]. 

Geometric numerical integration is concerned with devel- 
oping numerical integrators that preserve geometric features 
of a system, such as invariants, symmetry, and reversibility. 
In particular, variational integrators are geometric numerical 
integrators for Lagrangian or Hamiltonian systems, con- 
structed according to Hamilton's principle. They have de- 
sirable computational properties of preserving symplecticity 
and momentum maps, and they exhibit good energy behav- 
ior [11]. A variational integrator is developed for Lagrangian 
or Hamiltonian systems evolving on the two-sphere in [7]. It 
preserves both the underlying symplectic properties and the 
structures of the two-sphere concurrently. 

A variational integrator on for the controlled dynamics 
of a spherical pendulum can be written in a backward-time 
integration form as follows: 



qk 



2m/2 



2\ 1/2 



huJk+1 



2mP 



2mP 



2mP 



Ml 



k+l-, 



qu+i, (16) 



(17) 



where the constant > is time step, the subscript k denotes 
the value of a variable at the time tk = kh, and Mk = 
ml'^{-k^ujk - kqqd x qk). For given ((^/e+i, ^/c+i), we first 
compute Mfe+i. Then, q^ is obtained by ([16]), followed by 
M/e, and uJk is computed by ( [Tt] ). This yields an explicit, 
discrete inverse flow map T^^{{qk+i,uJk+i)) = {qk.^k)- 




Fig. 1. Stable manifold to (q,uj) = (—63,0) represented by 
{T~^{B§)}t>o for several values of t. One hundred points of B§ in 
the stable eigenspace to (—63,0) are chosen with 6 = 10~^, and they 
are integrated backward in time. Each trajectory is illustrated on a sphere, 
where the magnitude of angular velocity at each point is denoted by color 
shading (red: || 

^ max? 

blue: llc^llmin — 0)- 

3) Visualization: We choose 100 points on the surface of 
Bs with 6 = 10 and each point is integrated backward 
using ([16]), ( p7| ) with timestep h = 0.002. The resulting 
trajectories are illustrated in Fig. [T] for several values of 
t. Each colored curve on the sphere represents a trajectory 
on TS^, since at any point q on the curve, the direction of 
q = uj X q is tangent to the curve at q, and the magnitude of 
q is indirectly represented by color shading. 

We observe the following characteristics of the stable 
manifold VFs(— 63,0) of the inverted equilibrium: 

• The boundary of the stable manifold 63, 0) C TS^ 
parameterized by t is circular when projected onto S^. 

• Each trajectory in Ws{—es^O) is on a great circle, 
when projected onto S^. According to the closed loop 
dynamics ([T]), and the given initial condition at the 



surface of Bs, the direction of uj is always parallel to uj. 
Therefore, the direction of uj is fixed, and the resulting 
trajectory of g is on a great circle. This also corresponds 
to the fact that the eigenvalue Ai for the first mode 
representing the rotations about the first axis is equal 
to the eigenvalue A2 for the second mode representing 
the rotations about the second axis at ([T0|, i.e. the 
convergence rates of these two rotations are identical. 

• The angular velocity decreases to zero as the direction 
of the pendulum q converges to —63. 

• The stable manifold Ws {—es^O) may cover multiple 
times if t is sufficiently large, as illustrated at Fig. |l(f)| 
Therefore, at any point G S^, we can choose uj such 
that (q^uj) lies in the stable manifold W^{—es,0) (the 
corresponding value of uj is not unique, since if it is 
sufficiently large, q can traverse the sphere several times 
before converging to —63). This is similar to kicking 
a damped spherical pendulum carefully such that it 
converges to the inverted equilibrium. 

III. 3D Pendulum 

A 3D pendulum is a rigid body supported by a frictionless 
pivot acting under a gravitational potential. This is a gener- 
alization of a planar pendulum or a spherical pendulum, as 
it has three rotational degrees of freedom. It has been shown 
that a 3D pendulum may exhibit irregular maneuvers [12]. 

We choose a reference frame, and a body-fixed frame. 
The origin of the body-fixed frame is located at the pivot 
point. The attitude of a 3D pendulum is the orientation of 
the body-fixed frame with respect to the reference frame, and 
it is described by a rotation matrix representing the linear 
transformation from the body-fixed frame to the reference 
frame. The configuration manifold of a 3D pendulum is the 
special orthogonal group, SO (3) = {R e R^^'^lR^R = 
I,det[R] = 1}. 

The equations of motion for a 3D pendulum are given by 

JQ Q X JQ = mgp x R^es + u, 
R = R^^ 

where the matrix J G R^^^ is the inertia matrix of the 
pendulum about the pivot, and p G is the vector from 
the pivot to the center of mass of the pendulum represented 
in the body-fixed frame. The control moment at the pivot is 
denoted by G M^. 

A. Control System 

Several control systems have been developed on 
SO (3) [6], [4], [13]. Here, we summarize a control system 
to stabilize a 3D pendulum to a fixed desired attitude Rd G 
SO (3). Consider an attitude error function given by 

^{R,Rd) = ltT[{I-R^R)G], 

for a diagonal matrix G = diag[^i, ^2, ^3] ^ M^^^ with 
91^92^93 > 0. The derivative of this attitude error function 



with respect to R along the direction of SR 
is given by 

1 



Rfj for 7/ G 



Br^{R, Rd)-SR 
1 



{GR'^R-R' RdCy •r] = eR-r], 

where we use the property that tr[xA] = —x ■ {A — 
for any x eR^.Ae R^""^. The vee map, V : 5o(3) R^, 
denotes the inverse of the hat map. An attitude error vector 
is defined as cr = \{GR^R - R^RdG) e M^. For positive 
constants k^^kR, we choose the following control input: 

u = —kRCR — kf^ft — mgp x R^ e^,. 

The corresponding closed loop dynamics are given by 

m = -Vtx JVt- kRCR - knn, (18) 
R = Rn. (19) 

This system has four equilibria: in addition to the desired 
equilibrium (i?^,0), there exist three other equilibria at 
(i?(^ exp(7rei, 0), 0) for i G {1,2,3}, which correspond to 
the rotation of the desired attitude by 180° about each body- 
fixed axis. 

The existence of additional, undesirable equilibria is due to 
the nonlinear topological structure of SO (3), and it cannot be 
avoided by constructing a different control system (as long as 
it is continuous). It has been shown that it is not possible to 
design a continuous feedback control stabilizing an attitude 
globally on S0(3) [14], [15]. 

The stability of the desired equilibrium can be studied by 
using the following Lyapunov function, 

v = ^n-jn^kR^{R,Rd). 

In this paper, we analyze the stability of each equilibrium by 
linearizing the closed loop dynamics to study the equilibrium 
structures more explicitly. 

B. Linearization 

A variation in S0(3) can be expressed as [16]: 

R' = R exp(e7)) , = Q ^ eSQ, (20) 

for r],5Q G M^. The corresponding infinitesimal variation of 
R is given by SR = Rfj. Substituting this into (19), 

Rhfj -^Rfj = RfjCl + R6Q. 

Using the property xy — yx = x x y for any x, G 
can be rewritten as 

rj = SQ — CIt]. 
Similarly, by substituting ( [2Q| ) into ([18]), we obtain 

jdn = -sn X jn-nx jsn 

-\' 



^ this 



(21) 



-kR{GR^Rfi + fjR^RdG) - knSn, 
1 



where H = ir[R^ RdG]I - R^RdG G M^><^ and we used 
the property, xA + A^x = tT[A]I — A for any x eR^.A e 
]^3x3 Equations (21), ([22]) can be written in matrix form as 



-n I 

-\kRj-'^H J-\jh - nj - kni) 



Ax. 



(23) 



This corresponds to the linearized equation of motion of ([18]), 
([T9J. 

C. Equilibrium Solutions 

We choose the desired attitude as Rd = L In addition to 
the desired equilibrium (/, 0), there are three additional equi- 
libria, namely (exp(7rei), 0), (exp(7re2), 0), (exp(7re3), 0). 
We study the eigen- structure of each equilibrium using the 



linearized equation (23). We assume that 



J = diag[3,2,l]kgm^ G = diag[0.9, 1, 1.1], kR = kn = l. 

1) Equilibrium (/, 0).' The eigenvalues of the matrix A at 
the desired equilibrium (/, 0) are given by 

Ai,2 = -0.1667 ±0.5676z, 
A3,4 = -0.25±0.6614z, 
A5,6 = -0.5±0.8367i 

This equilibrium is an asymptotically stable focus. 

2) Equilibrium (exp(7rei), 0).- At this equilibrium, the 
eigenvalues and the eigenvectors of A are given by 



(24) 



Ai = 


-0.7813, 


Vi 


= ei 


- 0.781364, 


A2 = 


-0.5854, 


V2 


= 62 


- 0.585465, 


A3 = 


-1.0477, 


V3 


= 63 


- 1.047766, 


A4 = 


0.4480, 


V4 -- 


= eH 


h 0.448064, 


A5 = 


0.0854, 


V5 -- 


= e2H 


h 0.085465, 


A6 = 


0.0477, 


V6 -- 


= e3H 


h 0.047766. 



(22) 



Therefore, this equilibrium is a saddle point, where three 
modes are stable, and three modes are unstable. 

3) Equilibrium (exp(7re2), 0).' At this equilibrium, the 
eigenvalues and the eigenvectors of A are given by 

Ai = -0.3775, vi=ei- 0.377564, 

A2 = -1, ^2 = 62-65, (25) 
A3 = -0.9472, V3 = e3 - 0.947266, 
A4 = -0.0528, ^4 = 63 - 0.052866, 
A5 = 0.0442, ^5 = ei + 0.044264, 

A6 = 0.5, Vq = 565. 

Therefore, this equilibrium is a saddle point, where four 
modes are stable, and two modes are unstable. 



4) Equilibrium (exp(7re3), 0); At this equilibrium, the 
eigenvalues and the eigenvectors of A are given by 



Ai = 

A2 = 
A3 = 
A4 = 
A5 = 
A6 = 



-0.0613, 
-0.2721, 
-0.1382, 
-0.3618, 
-1.5954, 
0.5954, 



Vq = 



= ei 
= ei 

= 62 
= 62 
= 63 
62 - 



- 0.061364, 

- 0.272164, 

- 0.138265, 

- 0.361865, 

- 1.595466, 
0.595466. 



(26) 



Therefore, this equilibrium is a saddle point, where five 
modes are stable, and one mode is unstable. 

D. Stable Manifolds for the Saddle Points 

The eigen- structure analysis shows that there exist multi- 
dimensional stable manifolds for each saddle point. They 
have zero measure as the dimension of stable manifold is 
less than the dimension of TS0(3). But, the existence of 
these stable manifolds may have nontrivial effects on the 
attitude dynamics. 

We numerically characterize these stable manifolds using 



backward time integration, as discussed in Section II-D 



The stable eigenspace for each saddle point can be written 

as 

Ef,,(exp(^ei),0) = {{R,n) e TS0(3) | 
R = exp(7rei) exp(aiei + 0^262 + 0^363), 
n = -0.7813ai6i - 0.5854^262 - 1.0477^363 for G R] 

Ef,,(exp(^e2),0) = {{R,Q) G TS0(3) | 
R = exp(7re2) exp(aiei + 0^262 + {as + 0^4)63), 
Q = -0.37ai6i - 0^262 - (0.94Qf3 + 0.05^4)63 for G R 

EU^M^es)^0) = {{R,n) e TS0(3) | 
R = exp(7re3) exp((ai + ^2)61 + (0^3 + 0^4)62 + 61^563), 
n = -(0.06ai + 0.27^2)61 - (0.13^3 + 0.36^4)62 
- 1.590^563 for ai e M}, 
We define a distance on TS0(3) as follows: 

^TS0(3)((i^l, ^1), {R2.^2)) = ^^{Rl.R2) + 11^1 - ^2 ||. 

A variational integrator for the attitude dynamics of a rigid 
body on SO (3) is developed in [16], [17]. It can be rewritten 
in a backward integration form as follows: 

h 



^(H/c+i — -Mfe+i) = JdFk - Jd, 



life = Ffellfe+i 



Rk — Rk+lF]^ , 

h 

-FkMk+i 



(27) 
(28) 
(29) 



where — Uk-\-mgpxR-^ 63 G M is the external moment, 
Il/e = Jftk ^ R^ is the angular momentum. The matrix 
Jd G R^^^ denotes a non-standard inertia matrix given by 
Jd = ^tr[J]/ — J, and the rotation matrix Fj^ G SO (3) 
represent the relative attitude update between two integration 
time steps. For given 11^+1 ), we first compute M^+i, 

and solve ( [27] ) for F^. Then, Rk is obtained by ( [28] ), and 11^ 
is computed by ([29]). This yields a discrete inverse flow map. 




(a) Visualization on sphere (b) Rotation angle /3 (deg), and rota- 
tion rate /3 (rad/sec) 

Fig. 3. Visualization of an attitude maneuver: R(t) = exp(^(f)e3) for 
< t < 1, where ^(t) = |^(sin — 1). This maneuver corresponds to 
a rotation about the 63 axis by 30° to R(l) = I. The trajectory of the 
i-th column of R{t) representing the direction of the i-th body-fixed axis 
is illustrated on a sphere for i G {1, 2,3} (left). As the third body-fixed 
axis does not move during this maneuver, it is represented by a single point 
along the es axis on the sphere. The direction of R(t) is tangent to these 
curves, and the magnitude of R(t) is denoted by color shading, according 
to the magnitude of the rotation rate (right). 



1) Visualization (9/ VF5(exp(7rei), 0).' In [18], a method 
to visualize a function or a trajectory on SO (3) is proposed. 
Each column of a rotation matrix represents the direction 
of a body-fixed axis, and it evolves on S^. Therefore, a 
trajectory on SO (3) can be visualized by three curves on 
a sphere, representing the trajectory of three columns of a 
rotation matrix. The direction of the angular velocity should 
be chosen such that the corresponding time-derivative of the 
rotation matrix is tangent to the curve, and the magnitude 
of angular velocity can be illustrated by color shading. An 
example of visualizing a rotation about a single axis is 
illustrated in Fig. [3] 

We choose 112 points on the surface of Bs C 

10- 



£^f^^(exp(7rei), 0) with S = 10 ^, and each point is inte- 
grated backward using ([16]), (Vf) with timestep h = 0.002. 



The resulting trajectories are illustrated in Fig. |2] for several 
values of t. 

In each figure, three body-fixed axes of the desired attitude 
Rd = [ei, 62, 63], and three body-fixed axes of the additional 
equilibrium attitude exp(7rei) = [61,-62,-63] are shown. 
From these computational results, we observe the following 
characteristics on the stable manifold VF5(exp(7rei), 0): 

• When t < 15, the trajectories in Ws(exp(7rei), 0) 
are close to rotations about the third body-fixed axis 
63 to exp(7rei). This is consistent with the linearized 
dynamics, where the eigenvalue of the third mode, 
corresponding to the rotations about 63, has the fastest 



^d (^/e+l^n/e+l) {Rk.Uk). 



convergence rate, as seen in (24). 

When t > 15, the first mode representing the rotations 

about 61 starts to appear, followed by the second mode 

representing the rotation about 62. This corresponds to 

the fact that the first mode has a faster convergence rate 

than the second mode, i.e. |Ai| > IA2I. 

As t is increased further, the third body-fixed axis leaves 

the neighborhood of —63, and it exhibit the following 

pattern: 



O^^^^^^^^ ^^^^^^^^ fea ^ 

O O U 
T— 63 T— 63 T— 63 



(a) t = 11 (sec), \\Q\\ 
0.06 (rad/s) 



(b) f = 12 (sec), \\Q\\ 
0.17 (rad/s) 



(c) t = 13 (sec), \\Q\\ 
0.50 (rad/s) 



(d) t = 14 (sec), \\Q\\ 
1.42 (rad/s) 




(e) f = 15 (sec), 
3.93 (rad/s) 




(f) t = 16 (sec), ||f7||max 

10.67 (rad/s) 




(g) t = 17 (sec), 1 1 f7 1 1 max 

29.00 (rad/s) 




(h) t = 18 (sec), 
78.84 (rad/s) 



Fig. 2. Stable manifold to (exp(7rei), 0) = ([ei, —62, —63], 0) represented by *(B§)}t>o with S = 10 ^ for several values of t. 



• As t is increased further, nonlinear modes become 
dominant. The trajectories in Ws{ex.p{7re2)^0) almost 
cover SO (3). This suggests that for any initial attitude, 
we can choose several initial angular velocities such that 
the corresponding solutions converges to exp(7re2). 

3) Visualization of lVs(exp(7re3), 0); Similarly, we 
choose 976 points on the surface of Bs C £^f^^(exp(7re3), 0) 
with ^ = 10~^, and each point is integrated backward using 

0.002. The resulting trajectories 




• The stable manifold Wg {ex.p{7rei),0) covers a certain 
part of SO (3), when projected on to it. So, when 
an initial attitude is chosen such that its third body- 
fixed axis is sufficiently close to —63, there possibly 
exist multiple initial angular velocities such that the 
corresponding solution converges to exp(7rei) instead 
of the desired attitude Rd = L 

2) Visualization of iys(exp(7re2), 0).- We choose 544 
points on the surface of Bs C £^f^^(exp(7re2), 0) with 
S = 10~^, and each point is integrated backward using ( 16), 
^Tl) with timestep h = 0.002. The resulting trajectories are 
illustrated in Fig. |4] for several values of t. 

In each figure, three body-fixed axes of the desired attitude 
Rd = [ei, 62, 63], and three body-fixed axes of the additional 
equilibrium attitude exp(7re2) = [—61,62,-63] are shown. 
From these computational results, we observe the following 
characteristics on the stable manifold Ws(exp(7re2), 0): 

• When t < 12, the trajectories in VK5(exp(7re2), 0) is 
close to the rotations about the second body-fixed axis 
62. As t increases, rotations about 63 starts to appear. 
This corresponds to the linearized dynamics where the 
second mode representing rotations about 62 has the 
fastest convergence rate, followed by the third mode at 
(|25J. 



(T6|), ([17]) with timestep h 
are illustrated in Fig. |5] for several values of t. 

At each figure, three body-fixed axes of the desired attitude 
Rd = [61, 62, 63], and three body-fixed axes of the additional 
equilibrium attitude exp(7re3) = [-61,-62,63] are shown. 
From these computational results, we observe the following 
characteristics on the stable manifold VFs(exp(7re3), 0): 

• When t < S, the trajectories in Ws(exp(7re3), 0) are 
close to the rotations about the third body-fixed axis 
63. This corresponds to the linearized dynamics where 
the fifth mode representing rotations about 63 has the 



fastest convergence rate given in ([26]). 
• The rotations about 63 are still dominant, even as t is 
increased further. For the given simulation times, all 
trajectories in iy5(exp(7re3), 0) are close to rotations 
about 63. 

IV. Conclusions 

Stable manifolds of saddle points that arise in the closed- 
loop dynamics of two pendulum models are characterized 
numerically, and several properties are observed. Although 
the analytical and computational results have been presented 




0.03(rad/s) 0.09 (rad/s) 0.25 (rad/s) 0.69 (rad/s) 




(e) t = 15 (sec), ||17||max 

1.69 (rad/s) 




(f) t = 16 (sec), ||f7||max 

3.37 (rad/s) 




(g) t = 17 (sec), 1 1 17 1 1 max 

7.01 (rad/s) 




-63 



(h) t = 18 (sec), 1 1 17 1 1 max 

18.22 (rad/s) 



Fig. 4. Stable manifold to (exp(7re2), 0) = ([— ei, 62, — es], 0) represented by *(B§)}t>o with 5 = 10 ^ for several values of t. 




(a) t = 8 (sec), ||17|| 
0.224 (rad/s) 




(b) f = 9 (sec), \\n\\ 
1.09 (rad/s) 




(c) t = 10 (sec), \\n\\ 
4.26 (rad/s) 




(d) t = 14 (sec), \\n\\ 
222.99 (rad/s) 



Fig. 5. Stable manifold to (exp(7re3), 0) = ([— ei, —62, es], 0) represented by {T *(Bs)}t>o with 6 = 10 ^ for several values of t. 



for a spherical pendulum and a 3D pendulum, the methods 
presented naturally extend to any closed loop attitude control 
system with configurations in either or SO (3). 
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